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ABSTRACT 

We investigate the effect of coronal temperature on the formation process 
of solar chromospheric jets using two-dimensional magnetohydrodynamic simu¬ 
lations of the region from the upper convection zone to the lower corona. We 
develop a new radiative magnetohydro dynamic code for the dynamic modeling of 
the solar atmosphere, employing a LTE equation of state, optically thick radia¬ 
tive loss in the photosphere, optically thin radiative loss in the chromosphere and 
the corona, and thermal conduction along the magnetic field lines. Many chro¬ 
mospheric jets are produced in the simulations by shock waves passing through 
the transition region. We find that these jets are projected farther outward when 
the coronal temperature is lower (similar to that in coronal holes) and shorter 
when the coronal temperature is higher (similar to that in active regions). When 
the coronal temperature is high, the deceleration of the chromospheric jets is 
consistent with the model in which deceleration is determined by the periodic 
chromospheric shock waves. However, when the coronal temperature is low, the 
gravitational deceleration becomes more important and the chromospheric jets 
approach ballistic motion. 

Subject headings: magnetic fields — magnetohydrodynamics — Sun; atmosphere — 
Sun: chromosphere — Sun: transition region 
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Introduction 


Various chromospheric jets are observed in the solar atmosphere. In quiet regions 
and coronal holes, we observe classical (type I) spicules at the limb, which have maximum 


km/s (Beckers Il97 

2). Recentlv, the existence of more violent “tvpe H” spicules is reported 

(De Pontieu et al. 

2(1071 

). The tvpe I and H spicules are considered to have their on-disk 

counterparts of mottles ( 

Rouppe van der Voort et al. 2007) and the rapid blue-shifted events 

(RBEs) (De Pontieu et al. 2009: Rouppe van der Voort et al. 2009: McIntosh & De Pontieu 


20091) . respectively. The dynamic hbrils observed on the disk in active regions have 


maxi mum lengths of severa l Mm, lifetimes of 2-8 m in , and maximum velocities of 


Hansteen et al. 

2006; 

De Pontieu et al. 

2007a) 


Rouppe van der Voort et ah 


0-40 


2007h 


reported that the dynamic hbrils are considered to be the active region counterpart of the 
quiet-sun mottles. 


Theoretical s 


see 


Sterling 


2000 


udies have suggested various models for these chromospheric jets 
for a review). For th e dynamic hbrils, the periodic acoustic wave 


model is in agreemen t with observations flHansteen et al 


2006 


De Pontieu et al. 


Heggland et al. 


2007a : 


20071) . particularly in terms of the correlation between deceleration and 


maximum velocity. The origin of spicul es (types I an c 


candi dates, such as acoustic waves (e.g.. 




1969 


Hollweg et al. 


1982 


Heggland et al. 


Hollweg 


2009 


Takasao et al. 


I^udoh fc ShibatJ^ll999h . and magnetic reconnection (e.g.. 


1982 


II) is still under de 


Snematsn et al 


Date. Many 


19821) . Alfve n waves 


Uchida 


20131) have been suggested. Most of these models 


include a process by which the transition region is lifted by the shock wave propagating 


upward into t 


(iHollweg 


le corona. This process is called the shock-transition region interaction 


19821) . Several multi-dimensio nal simulations with the sophisticated mode 


convective motion have been reported (IHansteen et al. 


2006; 


Martmez-Svkora et al. 


ing o f 


2009 
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Heggland et al.ll201l[) . Although some of these studies is intended to be the representative 
for quiet regions, the simulated chromospheric jets was significantly shorter than the 
observed spicules. 


Using idealized one-dimensional hydrodynamic simulations, 


Shibata fc Suematsu 


(119821) explained why the spicules are long in the coronal hole. The key process is the 


amplification of the strength of the shock wave in the solar chromosphere. Based on 
their idea, we investigate the effect of the coronal temperature on the scale of the solar 
chromospheric jets using two-dimensional magnetohydrodynamic simulations with more 
realistic physical processes, to provide a unified perspective on the quiet region spicules, the 
coronal hole spicules, and the active region dynamic fibrils. 


2. Numerical Model 


We develop a new radiation magnetohydrodynamic code for the dynamical modeling 
of the solar atmosphere. The equations solved by this code include gravity, a LTE equation 
of state, the radiative cooling, and thermal conduction along the magnetic field line. 


The solar atmosphere is filled with strong shock waves with a wide range of plasma 
beta. A robust numerical scheme is required for the simulation of these regions. We employ 
a higher-order magnetp hydrodynamic scheme based on the constrained tra nsport inethod 


(Evans & Hawlev 

1988 

) with divergence-free magnetic field reconstruction ( 

Balsara 

2QQ9) 


The equations are solved in the conservative form. We do not consider explicit viscosity or 
magnetic diffusivity. 


We assume the local thermodynamic equilibrium for the equation o: 
most abundant elements are accounted with the metal mass fraction of 


(120061 ). The formation of molecular hydrogen is not considered. 


state. The six 


Asplund et al. 
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The radiative cooling term is evaluated as a combination of the two different 


approximations. In the photosphere and lower chromosphere, the radiative transfer 
equat ion is solved directly for the better modeling of the conve ctive motion flNordhmd 


19821 ). The short characteristic method (IKunasz fc Auer! Il988l) is employed with the 


A4 quadrature in 


Rosseland mean opacity 


opacity of 


Ferguson et ah 


Carlson! fll963h. The gray app roximation is assumed with the OPAL 


Iglesias fc Rogers 


1996!) in the high-temperature region and the 


(12005!) in the low-temperature region. In the upper chromosphere 


and corona, the optically thin radiative loss function is employed to estimate 


cooli ng. The loss function is calculated from the CHIA NTI atomic database flLandi et al 


20121) . with the extension bv !Goodman fc Judge! (120121) for the chromospheric temperature. 


he radiative 


These photospheric and chromospheric cooling terms are switched by the function of the 
column mass density. In addition to these terms, we include an artihcial heating term which 
becomes active only when the temperature falls below 2,500 K. 

Spitzer-type anisotropic thermal conduction is treated by second-order op erator 


splitting and the sec ond-order super-time 


stepp ing method (IMever et al. 


2012!) . The 


20071) is employed to maintain the monotonicity 


flux-limiting method (ISharma fc Hammett 
of the temperature. 

The two-dimensional numerical domain in the XZ-plane extends from 2 Mm below the 
surface to 14 Mm above. The horizontal extent is 18 Mm. The uniform grid size is 42 km 
in the horizontal direction and 32 km in the vertical direction. 

The bottom boundary condition is the “open” boundary. The total (gas plus magnetic) 
pressure equilibrium is assumed. The entropy of the upward plasma is fixed. The upward 


velocity held is slowly damped to the 
mass flux by the downward flow (e.g.. 


lorizon 

;allv 

Bercik 

2002) 


ally uniform upflow which complements the 


2002! ). This weak damping is not important 


in the present results but able to suppress the unrealistic instability which appears after 













































6 


the long time integration. The downward plasma evolves adiabatically and freely. The top 
boundary has a free-slip condition and is open for vertical flow. The density is extrapolated 
exponentially by the scale height at the boundary. Our two-dimensional simulation with a 


Leenaarts et ah 

2011; 

Heffdand et al. 

2011) 


20111 ). To prevent the corona from being cooled 


by radiation and thermal conduction, the conductive flux is introduced through the top 
boundary to heat the corona. In this study, the conductive flux is adjusted to preserve the 
different coronal temperature at the top boundary, T^. The magnetic field is assumed to be 
vertical at the top and bottom boundaries. The periodic boundary condition is employed 
horizontally. 


The initial condition is the plane-parallel atmosphere with a uniform vertical magnetic 
held of 3 G. The sufficiently relaxed convection is obtained after the integration of 5 solar 
hours with a temperature at the top boundary, T^, of 1 MK. We impose the uniform vertical 
magnetic held of 30 G on the relaxed atmosphere and integrate another 3 solar hours with 
the same Tc. Next, we change the temperature at the top boundary, Tc, to 2, 1, and 0.4 
MK to imitate the conditions of the active region, the quiet region, and the coronal hole, 
respectively. The simulations are integrated for another 1 solar hour after the change of Tc. 
We hnd that the coronal density and temperature reach statistical equilibrium within 20 
min. The data analyzed below is the last 30 min of each simulation. 


3. Results 

Figure [U shows snapshots from three simulations with the temperature at the top 
boundary, Tc, set to 2, 1, and 0.4 MK. We hnd the strong magnetic held concentration 
near X = 8 Mm. Because the magnetic held structure is essentially the same among 
these three simulations, we can concentrate on the ehect of the coronal temperature. The 









jet-like structures with various scales are formed in the chromosphere. These structures 
are basically driven by the shock waves, arising from convective motion, which pass 


through the 

transition region. T 

his process is consistent with both the one-dimensional 

simulations ( 

Ho 

weg et ah 

1982; 

Suematsu et ah 

1982 

) and multi-dimensional simulations 


flHeggland et af 


20111 ). We clearly hud that the simulation with a lower coronal temperature 


produces higher chromospheric jets. The density near the top of each jet and over the entire 
corona is smaller when is lower. The widths of the produced jets in these simulations are 
approximately independent of Tc if measured at a specihc height like Z = 3 Mm, but the 
jets become sharper near their tops in the lower simulations. 

We automatically detect the chromospheric jets in the simulati ons for the purpose o: 


quantitative comparison. The procedure is basically same as that in 


Heggland et ah 


11201 ih . 


The vertical motion of the transition region along the magnetic held line is htted by a 
parabola. The maximum length (i.e., the maximum height of transition region evaluated 
from the initial height), the lifetime, the maximum upward velocity, and the deceleration 
for each jet are calculated from the htting. The result of the statistical analysis is shown 
in Figure [21 We hnd a clear correlation between maximum length and maximum upward 
velocity for all three simulations (panels (a)-(c)). The lifetime roughly correlates with 
the maximum velocity for each simulation (panels (d)-(f)). The lifetime becomes longer 
and the correlation between lifetime and maximum velocity becomes stronger when the 
coronal temperature is lower. The correlation between deceleration and maximum velocity 
is relatively weak (panels (g)-(i)). The deceleration for Tc = 2 MK (panel (g)) is distributed 
over a relatively broad range, exceeding the gravitational deceleration on the solar surface. 
The deceleration for = 0.4 MK (panel (i)) is concentrated below the gravitational 
deceleration. 


The abovementioned dependence on the temperature at the top boundary, Tc, can be 



























explained as follows. When Tc is high (i.e., 2 MK) and the chromospheric jets are small, 
as in panels (a), (d), and (g) in Figure |21 the obtained correla t ions a re consistent with the 


‘shock deceleration hypothesis” suggested by 


Heggland et ah 


(120071) ■ According to this 


hypothesis, the deceleration is determined by the periodicity of the slow-moc 


e shock wave 


and can exceed the gravitational deceleration (see Eq. (1) in 


Heggland et al 


20071). The 


lifetime of the jet is determined by the period of the driving shock wave. The dotted lines 
in panels (a), (d), and (g) are the theoretical lines of parabolic motion with lifetimes of 
2, 3, and 5 min. These lines are well-consistent with the results from the simulation with 
Tc = 2 MK. Because the lifetime of chromospheric jet is equal to the period of driving 
acoustic wave in the “shock deceleration hypothesis”, the long lifetime in Tc = 0.4 MK 


short acoustic cutoff period in the lower solar atmosphere (e.g. 


3 not expectec 

under the 

Michalitsanos 

1973 

). The 


deceleration by the gravity becomes more important in this case. The resulting motion of 
taller chromospheric jets approaches free-fall. The dashed lines in panels (c), (f), and (i) 
indicate the theoretical lines of ballistic motion, which approximately explain the properties 
of simulated jets. Since the simulated jets in the case of Tc = 0.4 MK do not exactly follow 
ballistic paths, both of the gravitational and gas pressure gradient forces play a role in 
addition to the inclination of the magnetic held line. The case of Tc = 1 MK, shown in 
panels (b), (e), and (h), exhibits intermediate properties between these two extremes. 


To consider the abovementioned view on the determination of deceleration, we 
investigate the vertical force near the top of the chromospheric jets. The vertical component 
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of the equation of motion can be written as follows. 


DV., 


= D? + D!f + D^ + D 


G 


Dt 

Df = gobl 

= b- { -VP ) b 


D 




P 


( 1 ) 


D_i_= go- P|| = gob: 
p oz " 


Di = - 


1 

P 


J X B) 


Here, go is the gravitational acceleration at the solar surface, P is the gas pressure, B is 
the magnetic flux density, J is the current density, D / Dt = 8/dt + V • V is the Lagrangian 
derivative, and h = .6/ |.B| is a unit vector parallel to the magnetic held. If we assume 
that the transition region at the top of the chromospheric jets is the contact discontinuity, 
the deceleration of the jets can be approximated by —DVz/Dt near the transition region. 
The result is shown in Figure [3l Panels (a)-(c) show the deceleration produced by the 
total force parallel to the magnetic held versus the deceleration calculated by the 

motion of the chromospheric jets. Although there is some level of dispersion, the motion 
of the transition region is roughly explained by the forces parallel to the magnetic held. 
The possible candidates of the discrepancy are the ehect of the forces perpendicular to the 
magnetic held (P^, D^, and D^), the deviation of the transition region from the contact 
discontinuity by the thermal conduction and waves, and the htting and discretization errors. 
Panels (d)-(f) show the deceleration produced by the gravity parallel to the magnetic held 
Df. Df is more concentrated near the value of gravitational acceleration in the lower 
coronal temperature. Because the structure of the magnetic held is nearly independent 
of the coronal temperature, the temperature dependence in these panels are caused by 
the diherence of the height of produced jets. When the coronal temperature is higher, 
the maximum length of jets becomes small and the transition region moves near the root 
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of the magnetic field where the magnetic held is more inclined. Panels (g)-(i) show the 
deceleration produced by the gas pressure gradient force parallel to the magnetic held 
D^. In Tc = 0.4 MK case, the pressure gradient force always decreases the deceleration 
of the jets. However, in Tc = 2 and 1 MK cases, a part of the jets is decelerated by the 
pressure gradient force. This is consistent with the “shock deceleration hypothesis” where 
the deceleration is determined by the periodic shock wave. 


To understand why the scale of the chromospheric jets depends on the coronal 
temperature, we investigate the relationship between the coronal gas pressure and the 
maximum velocity of the produced jets as shown in Figure 01 The scale height of the gas 
pressure in the corona is sufficiently long and the gas pressure is almost constant across the 
transition region. Thus, we can use the coronal gas pressure, Pr, as t h e rep resentative value 


of the gas pressure near the transition region. 


Shibata fc Suematsul (119821 ) suggest that 


the density difference at the transition region causes the differences in the amplitude of 
the chromospheric shock wave and in the size of the chromospheric jets. For a linear wave, 
the energy flux of the slow magneto-acoustic wave is conserved during the propagation as 
SpV'^Cs = constant. Here S is the cross-section of the magnetic flux tube, Cg is the speed 
of the slow magneto-acoustic wave, and V is the velocity parallel to the magnetic held. 
This relation is reduced to H oc {SpCs)~^^'^ ~ The dependence on S does not 

appear in our current results because the magnetic held structure is almost the same in 
each of our three simulations. In fact, the shock wave is strong in the upper chromosphere, 
and the energy d issipation near the shock front violates the energy hux conservation. 


Ono et al. 


(11960) analytically derived the relation of H oc for a strong shock 


wave with the constant chr omospheric temperature and an adiabatic heat ratio of 5/3. 
Shibata fc Suematsul (119821 ) reported the relation of H oc in their one-dimensional 


adiabatic hydrodynamic simulation, which is very close to the analytic estimate for strong 
shock wave. We obtain the relation of K oc from the three simulations with different 
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coronal temperatures. This result indicates that another energy-loss process contributes 
to the formation of chromospheric jets in our simulations. The possible candidates are 
the radiative energy loss, the energy leakage by multidimensional effects, the thermal 
conduction, an d the latent heat of ionization. The shock-transition region interaction 


flHollweg 


19821) and the dynamic variation of the chromospheric temperature will also affect 


the result. We note that the amplitude of the maximum upward velocity depends only on 
the nonlinear propagation of acoustic waves in the photosphere and chromosphere and does 
not depend on the deceleration models discussed in the preceding paragraphs. 


Discussion 


We can compare our three simulations with different temperatures at the top boundary 
(Tc = 2,1, and 0.4 MK) to the observations of the active region dynamic hbrils, the quiet 
region spicules, and the coronal hole spicules, respectively. Panels (a), (d), and (g) in 
Figure [2] show behavior s imilar to the observation of the active region dynamic hbrils by 


De Pontieu et ah 


were investigated by 


(I2007al). The propertie s of the spicules in quiet 


Zhang et ah 


(120121) and 


Pereira et ah 


regions and coronal holes 


(120121) . Both studies reported 


that the spicules became longer in the coronal hole and shorter in the quiet region, which 
was consistent with our results. We note that the parabolic trajectory of jets in our study 
are consistent with the type I spicules but inconsistent with the type II spicules, which 


exhibit linear traiect ory. Because t 


observed wavelength (jPereira et al.l 


le top of the spicule was c onsiderably affected by the 


2014 : 


Skogsrud et ah 


20151) . further studies are required 


for the more quantitative comparison. We also note that = 0.4 MK is too low for 
coronal holes and = 2 MK is also too low for active regions. The regional difference of 
density and the spatial inhomogeneity of temperature and density will explain at least a 
part of this discrepancy. The magnetic held strength is also not the typical value for these 
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regions. However, if the chromospheric jets are driven by the acoustic waves from the lower 
atmosphere, we expect that the magnetic held does not affect the scale of the jets strongly. 


The numerical setting of the two-dimensional simulation by 


Heggland et ah 


fcniih 


was very similar to that in our study. The simulation with Tc = 1 MK is consistem 


Case B in their study. Their simulation employed the Bifrost code flGudiksen et al. 


with 


2niih 


with more sophistic ated treatment of radi a 


of our simulations. 


Martmez-Svkora et al. 


ion. This agreement supports the validity 


(120091) investigated the chromospheric jets in 


three-dimensional domain with the emerging magnetic hux. The typical length of their jets 
was approximately 1 Mm. The high coronal temperature produced by the energy release 
of the flux emergence, in addition to the rapid expansion of the magnetic held lines in the 
corona, is a possible explanation for why these jets were small. 


In this study, we neglect some physical processes which are important to the 


chromosphere, including non-equilibrium h ydrogen ionization (ICarlsson fc Stein 


more sophisticated treatments of radiation (iHayek et al. 


2010 


2002L 


Carlsson fc Leenaart? 


20121 ), 


and three-dimensionality. We believe that the qualitative picture is not affected. The 
dependence on the grid size is investigated and we hnd that the qualitative picture is not 
affected. 
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Fig. 1.— Snapshots from the three simulations with the different temperatures at the 
top boundary, Tc. Panels (a)-(c) show the temperature, and panels (d)-(f) represent the 
logarithm of mass density. The solid lines in panels (a)-(c) indicate the magnetic field 
lines. Panels (a) and (d) correspond to the simulation with T^. = 2 MK, panels (b) and 
(e) correspond to the simulation with Tc = 1 MK, and panels (c) and (f) correspond to 
the simulation with Tc = 0.4 MK. (An animation of this figure is available in the electronic 
edition of The Astrophysical Journal Letters.) 
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Fig. 2.— Statistical analysis of the produced chromospheric jets in the three simulations 
with different temperatures at the top boundary, Tc. Each cross sign corresponds to a chro¬ 
mospheric jet detected. See Section [3] for the description on the procedure for chromospheric 
jet detection. Panels (a), (d), and (g) correspond to Tc = 2 MK, panels (b), (e), and (h) 
correspond to Tc = 1 MK, and panels (c), (f), and (i) correspond to Tc = 0.4 MK. The three 
dotted lines in panels (a), (d), and (g) represent the theoretical lines of parabolic motion 
with lifetimes of 2, 3, and 5 min. The dashed lines in panels (c), (f), and (i) indicate the 
theoretical lines of ballistic motion. 
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Fig. 3.— Contribution of gas pressure gradient and gravitational forces on the deceleration 
of the produced chromospheric jets in the three simulations with different temperatures at 
the top boundary, Tc. Each cross sign corresponds to a chromospheric jet detected. Panels 
(a), (d), and (g) correspond to Tc = 2 MK, panels (b), (e), and (h) correspond to Tc = 1 
MK, and panels (c), (f), and (i) correspond to Tc = 0.4 MK. See Section[3]for the description 
on the dehnition of each variables. 
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Fig. 4.— Dependence of the maximum upward velocity of the chromospheric jets on the 
coronal gas pressure from the simulation of = 2 MK (shown as blue diamonds), = 1 
MK (shown as green crosses), and = 0.4 MK (shown as red triangle). The coronal gas 
pressure. Pc, is measured at the initial time of each chromospheric jet and averaged over the 
height range of 1-3 Mm above the transition region. The dashed line is the result of the 
least-square htting, 14 iax oc assuming a power-law relation. 



